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Abstract 

The equivalence of cooling to the gradient flow when the cooling step ric and the continuous flow step of 
gradient flow r are matched is generalized to gauge actions that include rectangular terms. By expanding 
the link variables up to subleading terms in perturbation theory, we relate Tie and r and show that the 
results for the topological charge become equivalent when rescaling r ~ tIc/( 3 — 15ci) where ci is the 
Symanzik coefficient multiplying the rectangular term. We, subsequently, apply cooling and the gradient 
flow using the Wilson, the Symanzik tree-level improved and the Iwasaki gauge actions to configurations 
produced with Nf = 2 -|- 1 -|-1 twisted mass fermions. We compute the topological charge, its distribution 
and the correlators between cooling and gradient flow at three values of the lattice spacing demonstrating 
that the perturbative rescaling r ~ nc/{3 — 15ci) leads to equivalent results. 
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1 Introduction 


Besides the interest by itself, the calculation of the topological properties of gauge held conhgurations is needed 
for several investigations in lattice QCD. These may involve a direct use of the topological charge in observables 
or its use as a measure of auto-correlations. The former, for example, includes the computation of the CP-odd 
form factor P 3 and subsequently the neutron electric dipole moment (nEDM) [1]. This would shed light on 
the question whether the value of the nEDM is zero or not and can therefore give hints of possible beyond 
the standard model physics. There is a number of smoothing techniques that could be applied to extract the 
topological charge Q, each one accompanied by its advantages and disadvantages mis]. The gluonic definition 
of the topological charge density in Euclidean time-space is given by 




( 1 ) 


with the gluonic field strength tensor and the totally antisymmetric tensor. The introduction of 

the gradient flow m 0 n] with its purturbatively proven renormalizability properties provides an attractive 
held-theoretic smoothing technique as compared to other techniques such as cooling and smearing, for which 
one can argue about the arbitrariness of their smoothing scale. The differential character of gradient-flow, 
however, makes it slower in comparison to other field-theoretic smoothers, such as cooling [7]. 

Recently it was demonstrated in Ref. [7] that using the Wilson action, gradient flow and cooling are 
equivalent if the gradient flow time r and the number of cooling steps ric are appropriately matched. By 
expanding the link matrices perturbatively in the lattice spacing a it was shown that at subleading order 
the two methods exhibit equivalence if one sets r = ric/S. This analytic result was verified by a numerical 
investigation of a number of observables such as the average action and the topological susceptibility confirming 
that the two procedures indeed produce equivalent results. This suggests that in cases where high statistics 
are needed such as, for example, for the evaluation of higher moments of the topological charge [5], instead of 
using the more expensive gradient flow, one can opt to employing cooling to evaluate quantities of interest. Of 
course in some applications, such as the scale setting through to-, where only a few hundreds of configurations 
are needed the computational cost is negligible and whether cooling or the gradient flow is used is not an 
important issue. 

Studies that utilize dynamical quark simulations such as those pursued by the European Twisted Mass 
Collaboration (ETMC) [SJlTnilll] make use of configurations produced with Symanzik improved gauge actions, 
such as the Iwasaki and the Symanzik tree-level improved actions |12j . It is interesting to extend the study of 
Ref. [7] to explore the use of Symanzik improved actions in the smoothing procedure. This choice will alter 
the relation between the scales r ~ ric/S since this depends on the choice of the smoothing action. We deliver 
the relation between gradient flow and cooling, by expanding the basic smoothing steps at subleading order in 
a for Symanzik improved actions. Subsequently we test the validity of the formula numerically using ETMC 
configurations produced with Nf = 2 -|- 1 -|- 1 twisted mass fermions and the Iwasaki gauge action. In addition 
to the Wilson action, we employ as smoothing actions the Symanzik tree-level improved and the Iwasaki 
actions enabling us to generalize the correspondence. We test the equivalence on the topological charge itself 
as well as on the average action and the susceptibility. We also examine the degree of correlation among the 
results obtained with cooling and the gradient flow through the correlation coefficient. All observables suggest 
that the two smoothers become equivalent after a few transient cooling steps. 

This article is organized as follows: In Section [21 we provide the relevant details regarding the production 
of the Nf = 2-I-1-I-1 configurations, in Section [31 we explain the different definitions of the topological 
charge density operators used for the calculation of the topological charge and in Section [H we provide a short 
description of the cooling and gradient-flow techniques for smoothing a gauge configuration in order to set the 
ground for their analytical comparison. We then compare the two smoothers by expanding the link variables 
perturbatively in a. In Section |S] we provide numerical evidence of this equivalence by evaluating a number of 
relevant observables. Finally, in Section [51 we summarize and conclude. 
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2 Configurations 


The gauge configurations are produced by the ETMC [5] using the Iwasaki improved action for the gluonic 
part 


= E {l-ReTr(C/i,"^iJ}+ci ^ {l - ReTr(t/ifJ} 


( 2 ) 




p,=^v 


with P = 2N/gQ, N = 3 and the plaquette and rectangular (1 x 2) Wilson loops. The Symanzik 

coefficients are set cq = 3.648 and ci = —0.331 and obey the relation cq + 8 ci = 1. The twisted mass fermion 
action at maximal twist is employed. The formulation provides automatic 0{a) improvement [131114) . infrared 
regularization of small eigenvalues and fast simulations with dynamical fermions. For the doublet of light 
quarks the action is 


S 


(i) 




= a^^X^'^\x)[Dw[U] +mo,i + , 


(3) 


where is the third Pauli matrix acting in the flavour space, tjiq i the bare untwisted light quark mass and 
g.1 the bare twisted light quark mass. The massless Wilson-Dirac operator is given by 


1 OT 


with the forward and backward covariant derivatives given by 
1 ^ 




U^{x)'tp{x + ajl) - ip{x) 


and V*il){x) = — 


Ul^{x — afl)ijj{x — ajjb) — il’ix) 


(4) 


(5) 


The fields “twisted basis” and are related to the fields in the physical basis through the 

transformations 


ip^''Hx) =-^{l + iT^j5)x^‘\x) and il}^''\x) = x^''\x)-^ {l + . 


( 6 ) 


Apart from the doublet of light quarks, we also include a twisted heavy mass-split doublet = (XcXs) for 
the strange and charm quarks. The associated action is expressed as 


S' 


[h) 


x^^\x^^\u = a'‘ E X^^\x){Dw[U] + TOo./i + g.5) , 


(7) 


with TOo,/i the bare untwisted quark mass for the heavy doublet, the bare twisted mass along the direction 
and the mass splitting in the direction. The heavy quark fields in the twisted basis are related to those 
in the physical basis through 

=-^ (a+ *t^75) ( 1 -f irVs) ■ (8) 

Unless stated otherwise, the quark fields will be understood as “physical fields”. The fermionic action in 
Eq. dS]) breaks parity and isospin at non-vanishing lattice spacing with the latter inducing a cut-off effect of 
0{a^) [14]. For more details on the twisted mass fermions see Ref. |9]. 

In order to test the equivalence between the two smoothing procedures we only need a single ensemble and 
a large number of configurations with a fine enough lattice spacing and relatively small pion mass. However, in 
order to investigate the behavior of observables as a function of the lattice spacing we include two additional 
ensembles, the pion mass of which is approximately the same as the one used for the more high statistics 
study. To this end, we selected the ensembles AGO.24, B55.32 and D45.32sc in the notation of Ref. m at 
three different lattice spacings so the continuum limit can be taken. The details of the ensembles can be found 
in Table (Tj 
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A60.24, /? = 1.90, a = 

= 0.094(1) fm, ro/a = 5.231(38) 

24^ X 48, L = 2.1 fm 

a/x 

0.0060 


No. of confs 

1160 


am-j^ 

0.17275(45)(23) 


LiJItt 

4.15 



0.362 GeV 

B55.32, /3 = 1.95, a = 

-- 0.082(1) fm, ro/a = 5.710(41) 

32^ X 64, L = 2.6 fm 


0.0055 


No. of confs 

4650 


am-j^ 

0.15518(21)(33) 


LuItt 

4.97 



0.372 GeV 

D45.32SC, /3 = 2.10, a 

= 0.064(1) fm, ro/a = 7.538(58) 

32^ X 64, L = 2.0 fm 

afj, 

0.0045 


No. of confs 

949 


am-j^ 

0.12087(40) 


LiJItt 

3.89 


rriT, 

0.368 GeV 


Table 1: Input parameters (/3, T, a/r) of our lattice calculation for the ensembles AGO.24, B55.32 and D45.32sc 
with the corresponding lattice spacing a, determined from the nucleon mass, and pion mass am,r in lattice 
units. 


3 Topological Charge 

3.1 Definition of the Topological Charge on the Lattice 

The topological charge of a gauge field is formally defined as the four-dimensional Euclidean integral over 
space-time 


Q = J d‘^xq{x) , (9) 

where the topological charge density q{x) is defined in Eq. ([T]). 

In practice, any valid lattice discretization of q{x) —>■ qL^x) leading to the right continuum expression 

of Eq. dH) can be used for the evaluation of the lattice equivalence of Eq. (|9|), given by 

Q = a^'^qLix). ( 10 ) 

X 

However, depending on the discretization of the operator qL{x) lattice artifacts affecting the total topological 
charge Q vary. Hence, we do not expect to obtain an exact integei0 value for the topological charge. Nev¬ 
ertheless, we expect that the total topological charge, for some definitions for the topological charge density, 
converge faster and are closer to an integer than that obtained by other definitions. To investigate the dif¬ 
ferent definitions we use a number of lattice discretizations. The simplest lattice discretization, which can be 
constructed is based on the simple plaquette, depicted pictorially in Fig. [T] 

= Im [U^j,{x)Ui,{x + ajl)Ul{x + av)Ul{x)] , (11) 

^Of course one can obtain an exact integer when applying the Atiyah-Singer index theorem |16II17| Q = n_ —n+ and employing 
the number of Dirac zero modes n± with positive (+) and negative (-) chiralities obtained with the Overlap-Dirac operator m- 
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with 


( 12 ) 


This is a computationally cheap definition which, however, leads to lattice artifacts of order 0{a^). Nev¬ 
ertheless, this is still an adequate definition having been used in several determinations of the topological 
susceptibility in the past [TH [20] . 

Indubitably, the most common definition of the topological charge density is the clover definition given by 

{G^'rC^'r} , (13) 

with G(('°''(a;) the usual clover leaf (second picture in Fig. |T|) defined as 


G/'r(^) = 

Uf^{x)Uu{x + afL)Ul{x ai')Ul{x) 


+ 

U^{x)Uj^{x — afL + av)Ul{x — afL)U^{x — afi) 


+ 

U^{x — afL)Ul{x — afi — av)Ufj,{x — afi — av)Uv{x — av) 


+ 

Ul(x — av)Ufj,{x — av)Uv{x + afi - av)Uj^{x) . 

(14) 


However, this definition still carries a leading correction term of G(a^). Hence, an improved definition of 
the topological charge density, which removes tree-level discretization errors and converges as 0(a'^) in the 
continuum limit is also considered. Such a definition, given in Refs. [DllTlIll], is 


= Coqf°^{x) + Cxq^£^\x) , 


(15) 


where qf^°'^(x) is the ordinary clover topological charge density in Eq. dT51) and q^£^^(x) is the clover-like 
operators where instead of squares we make use of horizontally-and vertically-oriented rectangular Wilson 
loops of size 2x1 and 1x2 respectively 

qT\^) = ^ WTr {G-*G^r } , ( 16 ) 

with 


G 


rect 


{x) 


— Ufi{x)Uu{x + ajl)Uv{x -\- a[l + av)U^^{x + 2av)Ul{x + av)Ul{x) 

Ui,{x)Uu{x + av)Uj^{x — ajl + 2av)Ul{x — ajl-\- a{')Ul{x — ajl)U^{x — afi) 

Uj^{x — afi)Ul{x — afi — av)Ul{x - aft - 2av)U^{x — afi — 2av)Ui,{x — 2ai>)Uu{x 
Ul{x — av)Ul{x — 2av)U^{x — 2aC')U,y{x + afi — 2ai>)Ui,{x + af — ai')U^{x) 
U^{x)U^{x + afi)U^{x + 2afi)U^{x + ai) + afi)U^{x + ai))Ul{x) 

Ui,{x)U^{x — aft + ai))t/l(x — 2afL + av)Ul{x — 2afL)U^{x — 2afL)U^{x — afi) 
Uj^fx — afL)Uj^{x — 2afi)Ul{x — 2afi — av)U^{x — 2af. — av)U^{x -aft- av)U,^{x 
Ul{x — av)U^{x — aC')Ufj,{x — ai> + afi)Uu{x + 2afi — av)Uj^{x + afi)U^{x) . 


av) 


av) 


(17) 


In order to remove the discretization error at tree-level one should use the Symanzik tree-level coefficients 
Cl = —1/12 and cq = 5/3. A diagrammatic representation of the three definitions of Gf^{x) (r =plaq, clov, 
rect) used in our investigation is provided in Fig. [T| 

Ultraviolet fluctuations of the gauge fields entering in the definition of e.g. the topological charge lead to 
non-integer values. Thus, methods to suppress these ultraviolet fluctuations are employed. Such techniques 
include cooling and the more recently introduced gradient flow. We examine both these techniques using, 
beyond the Wilson, the Symanzik tree-level improved and Iwasaki actions. 
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Figure 1: From left to right, we represent pictorially the plaquette operator used for the definition of the 
the ordinary clover 0^"' and the rectangle clovers ^ g^^h that ^ 


4 Equivalence of Cooling with Gradient-Flow 

We smooth out the ultraviolet fluctuations using the action in Eq. ©. The Symanzik coefficients must satisfy 
Co + 8ci = 1 and aside from this requirement, the value of ci can be chosen arbitrarily. The case of ci = 0 
corresponds to the ordinary Wilson action. In addition to the Iwasaki action we also consider the Symanzik 
tree-level improved action with ci = —1/12. Any discrepencies resulting from different smoothing actions are 
interpeted as lattice artifacts and are expected to vanish in the continuum limit. 

Smoothing a gauge link C/^(x) can be accomplished by its replacement by some other link that minimizes 
the local action. To this purpose it makes more sense to rewrite the gauge action of Eq. @ as 

Sg = ^ReTr{Ajj(a;) [/^(x)}-I-{terms independent of U^{x)} , (18) 

where X^x) is the sum of all the path ordered products of link matrices, called the “staples”, which interact 
with the link U^{x). The main components in the Wilson action are the plaquettes and thus the staples 
resulting from the square component of the action extend over 1x1 squares (in lattice units). For the 
rectangular part of the action the staples extend over rectangles of sizes 1x2 and 2x1. We can, therefore, 
write Xfj_{x) as 

XTx) = coXf:'Tx) + ciXT\x) , (19) 


with 


xT^ix) = 


E 


Uv{x)UTx + ai))Ul{x + ajl) + Ul{x — av)UTx — av)Uy(x — av + ajl) 


( 20 ) 
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and 


^ \^i,{x)Uv{x + ai')u^{x+ 2ai>)ul{x + av + aii)ul{x + afl) + 

+ Ul{x — av)Ul{x — 2av)U^{x — 2av)U,^{x — 2ai> + ajl)Uv{x — av + afi) 

+ ^ \^y{x)Ufj_{x + av)Ufj,{x + av + ajl)Ul{x+ 2ajl)Ul^{x + afl) + 

+ U},{x — av)U^{x — av)U^{x — av + ail)Ui,{x — av + 2afi)U^{x + ajl) 

+ \^‘l^(x — ajl)U,^{x — ajl)U^{x — ajl + av)U^{x + av)Ul{x + ajl) + 


+ ~ ajl)Ul{x — av — ajl)U^{x — av — ajl)U^(x — av)Uu{x — av + ajl) 


( 21 ) 


According to the above two equations, for a given link U^{x), the total number of plaquette and rectangular 
staples interacting with it is 6 and 18 respectively. 


4.1 Cooling 


Cooling is applied to a link variable 17^(x) G SU{N) by updating it, from an old value to Ujj™{x), 

according to the probability density 


P{U) oc exp 


lim /3^ReTrA^'l'(a;)t/^(a;) 
3 —¥(X> iV 


( 22 ) 


The basic step of the cooling algorithm is to replace the given link Ujj^'^jx) by an SU{N) group element, which 
minimizes locally the action, while all the other links remain unaltered. This is done by choosing a matrix 
U^™{x) G SU(N) that maximizes 

ReTr{[/“-(x)A^t(3,)}. (23) 

In the case of an SU(2) gauge theory, the maximization is achieved by 


TT-new 


(x) 


y^detXfj,{x) 


(24) 


For SU{N) the maximization can be implemented by using the Cabibbo-Marinari algorithm [23) : one has to 
iterate the maximization over all the SU{2) subgroups embedded into SU{N). 

We iterate this procedure so that all the links on all sites are updated. Such a sweep over the whole lattice 
is called a cooling step and will denote by ric the number of cooling steps performed. During the sweep the link 
variables, which have already been updated, are subsequently used for the update of the links still retaining 
their old value. 


4.2 Gradient flow 

The gradient flow is defined as the solution of the evolution equations HISJIS] 

Vf,{x,T) = -gl[d^^f,SGiV{T))]Vf,{x,T) 

V^{x,Q) = U^ix), 

where r is the total gradient flow time. In the above equation the link derivative is defined as 


d=,,^SG{U) 


a 

iY.T%iSG{U), 

a 


s=0 


(25) 


(26) 
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with 


Y^y.y) 


T“ if (y, i/) = (x,/i) 

0 if (y, i^) ^ (x, y), 


(27) 


and T°‘ (a = 1, • • • , — 1) the Hermitian generators of the SU{N) group. If we now set 

we obtain 


9ld.,^SG{U) = O],) - ^Tr - fit) . (28) 

The last equation provides all we need in order to smooth the gauge fields according to the Eqs. (123. Evolving 
the gauge fields via gradient flow requires the numerical integration of Eqs. (l25l) . This is performed using the 
third order Runge-Kutta scheme as explained in Ref. [B]. For the exponentiation of the Lie-algebra fields 
required for the integration, we apply the algorithm described in Ref. [24] . We investigate how the elementary 
integration step e affects our results and find that e = 0.02 is a safe option as this was also pointed out in 
Ref. [7]; we observe that smaller elementary integration steps give the same results. We therefore set e = 0.02 
for the integration step. 

4.3 Perturbative relation between cooling and the gradient flow 

Both cooling and gradient flow can be used to remove the ultraviolet fluctuations. Both should lead to the 
same topological properties provided that we are close enough to the continuum limit. Assuming that we are 
in the perturbative regime we can carry out a perturbative comparison in order to obtain an analytic relation 
between the scales involved in the two procedures following Ref. [7] where the relation 

r~nc/3, (29) 

was derived for the Wilson action. In this work we derive a more general expression of the form r = ric x /(ci) 
for smoothing actions that, in addition to the plaquette, also include a rectangular term. 

In the perturbative regime the link variables can be expanded as 

U^{x)^\ + iY,ul{x)T\ (30) 

a 


with uj)(a;) € R is assumed to be inhnitesimal. 

Using Eqs. (EOl) and (ED the plaquette and rectangular staples are written as 

XP'"'i(x):x6-l-*^<(x)r“ and X),-‘(x) 18 + * ^ z;“(x)r“ , (31) 

a a 

where w“{x) and v“{x) are infinitesimal quantities. The leading coefficients with values 6 and 18 appearing in 
the above equations are just the number of plaquettes and rectangles interacting with the link on which the 
gradient flow evolution is applied. We can, therefore, write the sum of staples (Eq. (flOll l as 

Xij_{x) ~ 6co -b 18ci -b ico ^ w“(a;)T“ -b ici ^ f“(a;)T“ , (32) 


and, subsequently, U^(a;) as 

Q^{x) ~ 6co -b 18 ci -b i ^ [(6co -b 18 ci)u“ (x) - (cow“(x) -b ciu“(x))] T“ . (33) 

a 

Hence, Eq. (1751) becomes 

9ldx,tiSG{U) = i ^ [(6co -b 1801)14“ (x) - (coic“(x) -b cii;“(x))] r“ . 

a 
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(34) 


Using the above expression, the evolution of the gradient flow can be approximated as 


(x, r + e) ~ r) - e [(6co + 18ci)u“(x, r) - (cow“(a;, r) + ci?;“(a;, r))] . (35) 

For the cooling procedure, one needs to consider that the link Ufj,{x) is substituted with the projection of 
over the gauge group. Namely, for the case of an SU{2) gauge theory this projection is manifested by Eq. (l24ll 
where we substitute (x) by Eq. (15^ . In the perturbative approximation this leads to H 


'(x) ~ 1 + i 


(cowg(x) +CiX°(x)) 
6 co + 18ci 


rpa 


The above update corresponds to the substitution 

(cow“(x) +ciu“(x)) 




6 co + 18ci 


(38) 


(39) 


Comparing Eqs. [35] and [39| we observe that the gradient flow would evolve the same as cooling if one chooses 
a step of e = l/(6co + 18ci). In addition, during a whole cooling step the link variables, which have already 
been updated are subsequently used for the update of the remaining links that await update; this corresponds 
to a speed-up of a factor of two. Therefore, the predicted perturbative relation between the flow time r and 
the number of cooling steps Uc so that both smoothers have the same effect on the gauge field is 


^ ric _ Uc 

3co + 9ci 3 — I5ci 

The cooling/gradient flow rescaling factors for our choice of actions are given in Table [2| 


(40) 


Smoothing action 

Co 

Cl 

Uc/r 

Wilson 

1 

0 

3 

Symanzik tr.level 

5 

3 

1 

12 

4.25 

Iwasaki 

3.648 

-0.331 

7.965 


Table 2: Leading order perturbative rescaling between the number of cooling steps and gradient flow time 
such that the two smoothing techniques are equivalent. These numbers are according to Eq. (|40|). 


An important question, which needs to be answered is how one tunes the smoothing parameters as the 
continuum limit is approached; this has been extensively discussed in Ref. [7] and we will briefly comment 
on how this is modified here. In practice, by applying the smoothing procedure on some configurations the 
ultraviolet (UV) properties of the theory up to some length scale Ag are modified by suppressing the UV 
fluctuations at smaller length scales. Eor this to be a viable procedure we need to show that by altering the 
UV part of the theory the continuum results remain unchangeable and, thus, the underlying physics does not 
depend on Ag. Thus, one needs to choose the length scale As, which for procedures like cooling often was 
taken arbitrarily; in other words the choice of the smoothing parameters such as ric in the case of cooling but 

^This can be derived easily for SU{2) where one can explicitly expand Eq. II24II . Making use of the Mercator series expansion 
of the logarithm we write 


detX^ (x) = (6co + 18ci)^det 

\ a 

Thus, the expansion of Eq. isi gives 
Xu 


(cpm^fa;) -hcmg fa:)) ^ 
6co + 18ci 


yMetXjTfx) (6co + 


^6co -I- 18ci -h ICO ^ w°{x)T‘" + ici ^ v°{x)T°‘'^ = 1 + i ^ 


T“ = (6co -h 18ci)2 {1 + 0 (Tr {r“}) -f O (a^)) . (36) 


{cow'^jx) + civ^jx)) 
6co + 18ci 
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also for other smearing techniques such as APE [^[55], HYP miiiH] and Stout PU [52] is not entirely clear. 
The gradient flow, on the other hand, provides a smoothing procedure where this length scale is quantified 
as discussed below. Namely, it has been shown that one can simply renormalize composite operators at fixed 
physical flow time with 


Xs-Vm, (41) 

with t = a^T being the gradient flow time in physical units. We can, therefore, translate the length scale Xs 
as a function of the cooling step ric according to the formula 


Xs ^ a 


Sric 


3 - 15ci 


(42) 


Given that the validity of Eq. (gni) is confirmed numerically, we end up with an expression for an associate 
length scale Xs for the case of cooling as well. One can also generalize this correspondance for the cases of 
other smoothers, such as the APE and stout smearing [30] . 

As an example we consider the continuum limit of the topological susceptibility which is used in this 
work. According to Refs. HE] one reads the topological susceptibility at a fixed value (in physical units) of 
As = = 0(0.Ifm) such that As is not too small so that discretisation effects are suppressed, as well as 

not too large so that the topological content of the gauge field is preserved. Practically, As should correspond 
to a plateau for the topological susceptibility which should be scale invariant. Hence, we extract the value 
of the topological susceptibility at fixed As for a sequence of lattice spacings and then extrapolate it in the 
continuum limit. 


5 Numerical Results 

5.1 Topological Charge 

We apply cooling and gradient flow on TV/ = 2+1 + 1 twisted mass fermions gauge configurations with /3 = 1.90, 
/3 = 1.95 and /3 = 2.10 using the Wilson (Eq. ([2]) with ci =0), Symanzik tree-level improved (Eq. ([2]) with 
Cl = —1/12) and Iwasaki (Eq. ([5]) with Ci = —0.331) actions. We measure the average action, as well as 
the plaquette (Eq. (fT51) '). clover (Eq. (IT51) 1 and improved (Eq. (HU)) definitions of the topological charge for 
every cooling step ric- Gradient flow is costlier and, thus, we take measurements for every Ar = 0.1 in units 
of gradient flow time (which corresponds to five integration steps for e = 0 . 02 ) instead of every integration 
step. We cover in total 60 — 100 cooling steps while for the gradient flow we fix the maximum gradient flow 
time according to the perturbative expression of Eq. (1401) and the maximum number of cooling steps. The 
cooling/gradient flow rescaling factors used are taken from Table [5] 

The behavior of the topological charge Q for single configurations as a function of ric and r is investigated 
for cooling and gradient flow, respectively for a given smoothing action and lattice spacing. In Eig. |2| we 
present the improved definition of the topological charge as a function of ric and t for four randomly chosen 
configurations. We show results for (3 = 1.90 and /3 = 2.10. For j3 = 1.90 we observe that the topological 
charge for a given configuration whithin the whole range of ric j (3 — 15ci) x r yield different values for cooling 
and gradient flow. The difference in the value of the topological charge is not surprising since the different 
smoothers have different lattice artifacts and do not need to agree at non-zero values of the lattice spacing. 
For j3 = 2.10 the values become closer as expected. Thus, as one approaches the continuum limit the two 
different procedures converge. We note that the topological charge itself is not the main quantity of interest. 
It provides only a measure on the fluctuations and an input for the topological susceptibility, which is the 
physically relevant quantity. In the next section we will thus focus on the relevant physical observables. In this 
section, we restrict the presentation to the topological charge. Another observation from the results shown in 
Fig. |2|is that for the Wilson and Symanzik tree-level improved actions the topological charge Q as a function 
of ric or T is not really constant. As can be seen in the left and middle panels of Fig. H the topological 
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ric or 3 X T Tie or 4.25 X r ric or 7.965 X r 

Figure 2: The improved definition of the topological charge as a function of ric for cooling and r rescaled 
by a factor of 3, 4.25 and 7.965 for gradient flow extracted with Wilson (left), Symanzik tree-level improved 
(middle) and Iwasaki (right) smoothing actions respectively. The different colors correspond to the four 
different configurations chosen randomly while filled and open symbols correspond to cooling and gradient 
flow respectively. Upper panel is for fj = 1.90 and lower panel for /3 = 2.10. 


charge obtains different values with increasing ric and r. This behavior, although still present, appears to 
be supressed for our finest lattices with fi = 2.10. Using the Iwasaki action, we observe that the topological 
charge fluctuates for ric G [0, 20 — 30] and then becomes completely stable no matter what the lattice spacing is. 
These results have been observed when applying cooling in previous studies and they comply with theoretical 
expectations from an, admittedly, semi-classical picture. Namely, at finite lattice spacing the lattice action 
deviates from its continuum limit with deviations that increase as the gauge fields become larger. Instantons 
have a scale parameter A, which enters non-trivially the action. As one decreases A, the gauge fields are 
expected to become larger modifying the gauge action as well. The lattice action can be written [SIEI] (on 
dimensional grounds) as 

5Lat(a, A) = 5cont {l + {a/\f a2 + {a/xf a^ + O (a/A)®} (43) 

with 02 = —1/5 for Wilson, 02 = 0, 04 = —17/210 for Symanzik and 02 = -1-2.972/5 for Iwasaki. Stable instan- 
ton solutions require a lattice action which increases by decreasing the scale parameter A. This requirement is 
fullfield only for the Iwasaki action and that is the reason why one observes stable topological charge. On the 
contrary, for the Wilson and Symanzik tree-level improved actions, the solutions are not stable; this is reflected 
in the fact that the values of the topological charge jump to different values. Nevertheless, stability sets in 
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as a —?► 0; this is visible for the case of /3 = 2.10 in Fig. [5] where we observe less changes in the value of Q. 



Figure 3: An example of the behavior of the topological charge for a single conhguration as a function of 
ric and r rescaled by 7.965 for cooling (open symbols) and gradient flow (filled symbols) for two different 
configurations. With the red circles we present the improved, with green diamonds the clover and with blue 
squares the plaquette definition of the topological charge. The smoothing has been performed with the Iwasaki 
action. 

Comparing results for the three different definitions of the topological charge density, we observe that for the 
improved case the topological charge converges closer and faster to a near integer value compared to the other 
two definitions. All three definitions for the three ensembles give topological charges, which converge to the 
same near integer as a function of the relevant smoothing scale. These two observations suggest that indeed 
the three topological charge definitions differ only due to lattice artifacts. Such a comparison is meaningful 
only if the topological charge acquires stability and hence, we consider the Iwasaki action. In Fig. [31 we 
observe that for the clover as well as for the improved definition, the topological charge converges faster than 
when the plaquette definition is used in particular in the case when cooling is performed. 

In Fig. HI we present an example of the time history (first 900 configurations) of the topological charge Q 
for gauge configurations that have been cooled using Uc = 50. We also include the time history when using 
the gradient flow for a step of t = nc/(3 — 15ci). Results are shown for (3 = 1.95 and /3 = 2.10 for the three 
gauge actions. As can be seen, the topological charge does not suffer from large autocorrelations and the 
time histories between cooling and the gradient flow are very similar. This similarity can be quantified by the 
calculation of the linear correlation coefficient, which is the topic of Section (AH 

Additionally, in Fig. lU we provide the histogram of the topological charge for both cooling and gradient 
flow for the three actions. We observe that the histograms exhibit nearly gaussian distributions in particular 
for the P = 1.95 ensemble where a large number of conhgurations is analyzed. As expected, the distributions 
using cooling or the gradient flow look very similar for all three actions and the associated gaussian fits fall 
on top of each other. This already points to the equivalence anticipated for the topological susceptibility. 

5.2 Average Action Density 

As a common scale for the two smoothing techniques we can use the action, the minimization of which defines 
both smoothers. Instead of looking at the action we consider the dimensionless average action density (^g) 
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Figure 4: The time history of the topological charge which has been extracted by cooling (blue dashed line) 
and gradient flow (red solid line) at Uc = 50 and the corresponding flow time for each different choice of 
smoothing action. In the upper row we present results for (3 = 1.95 and in the lower results for jS = 2.10. 
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In Fig. ini we present the average action density for /3 = 1.95 as a function of Uc and the perturbatively 
determined values of the gradient flow time, namely 3x, 4.25x, and 7.965 x r for the Wilson, the Symanzik 
tree-level improved and the Iwasaki actions, respectively. As expected from the findings of Ref. [7], for the 
Wilson action, the rescaling Uc = 3r leads to equivalent results for this quantity between gradient flow and 
cooling for small values of Uc and r. For instance for /? = 1.95 where our results are more accurate we find that 
for ric > 20 the average action for both procedures becomes the same. Our results show that a similar behavior 
is observed also for the other two actions. Namely, the average action density deviates for small values of the 
smoothing scales but for nc^30, for the Symanzik tree-level improved, and nc~50, for the Iwasaki action, they 
become equal. Similar behavior is also observed for /3 = 1.90 and [3 = 2.10 showing the equivalence of the 
two procedures in evaluating the average action density. Following Ref. [7] we define r(nc) as the gradient 
flow time r for which the average action density changes by the same amount as when ric cooling steps are 
performed. This function is evaluated by interpolating between the discrete gradient flow time steps with 
cubic splines. In Fig. [7] we report the function r(nc) for the three different actions, for our three different 
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Figure 5: The distribution (first row) of the topological charge for j3 = 1.95 and the accosiated gaussian fit 
(second row). In blue we present the distribution obtained via cooling at ric = 50 and in red the distribution 
obtained via gradient flow at r = 16.7, r = 11.8 and r = 6.3 for Wilson, Symanzik tree-level improved and 
Iwasaki actions, respectively. 


ensembles. We observe that for each action used the results are in agreement for the three ensembles giving 
the first indication that the equivalence between the gradient flow and cooling has a well-defined continuum 
limit. In addition to the functions T(nc) we also plot the lines r = nc/3, r = nc/4.25 and r = nc/7.965 for 
the Wilson, Symanzik tree-level and Iwasaki actions, respectively. Obviously these linear functions provide 
good approximations of r(nc) for each choice of action even for the ranges of ric where equivalence in Fig. [6] 
does not hold. Since the average action plays the role of a common scale between the two procedures and 
T(nc) has such a good agreement with the perturbative lines, there is little doubt that the approximation 
T(nc) = nc/(3 — I5ci) provides an adequate rescaling between ric and t with finite lattice spacing corrections 
playing an insignificant role. 


5.3 Topological Susceptibility 

In this section we examine results on the topological susceptibility defined as 

^ aW ■ 


(45) 


The topological susceptibility has been investigated extensively using several techniques such as smearing and 
cooling [221 EH [32] and recently determinations of y make use of the gradient flow [33] as well as the spectral 
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Figure 6: The average action density (Sg) as a function of the cooling step ric and the corresponding gradient 
flow time nc/(3 — 15ci) for /3 = 1.95 and the three smoothing actions Wilson, Symanzik tree-level (tr.l) and 
Iwasaki. 



Figure 7: The behavior of T(nc) as a function of ric for Wilson, Symanzik tree-level improved and Iwasaki 
smoothing actions. The lines corresponds to r = ndS, t = nc/4.25 and r = nc/7.965 for Wilson, Symanzik 
tree-level improved and Iwasaki actions, respectively. 


projectors method [331133 ■ The question we would like to address here is not the detailed determination of 
the topological susceptibility, which will be the subject of another followup paper, but rather its use as a 
comparison between cooling and gradient flow for the three actions considered in the previous sections. In 
Fig. |8]we show as a function of the number of cooling steps and the gradient flow time rescaled by the 

corresponding perturbative factor for /3 = 1.95. We do so for the three different actions used in the cooling 
and gradient flow procedure, namely the Wilson, the Symanzik tree-level improved and the Iwasaki actions, 
for the three lattice definitions of the topological charge density; to reveal the associated correspondence we 
collect the results for both procedures in the same plot. We observe that for a given action after a few cooling 
steps Uc < 10 or the equivalent gradient flow time r = nd d ~ 15ci) the susceptibilities computed using the 
plaquette or the clover definition of the topological charge density are almost indistinguishable. In general, such 
an agreement is not expected at finite lattice spacings and one might see deviations for very large statistics. 
However, for our current statistical accuracy both definitions give the same results and we thus only considered 
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Figure 8: The topological susceptibility units of tq computed using the three different definitions of the 
topological charge density, namely the plaquette, the clover and the improved definition, as a function of the 
cooling step and the associated gradient flow time. From left to right we show results for the Wilson, Symanzik 
tree-level improved and Iwasaki actions. The results when using the plaquette definition coincide with those 
obtained using the clover definition. 


the susceptibility based on the clover definition of the topological charge. The results in Fig. [8] also show the 
very good agreement between cooling and the gradient flow for the topological susceptibilities obtained using 
the same definition for the topological charge density and the same action. As a matter of fact for even a 
very small number of cooling steps i.e. ric ~ 5 and the corresponding gradient flow time r ~ 5/(3 — 15ci) 
the two values of the topological susceptibilities become the same. For larger number of cooling steps and 
the associated gradient flow times the two topological susceptibilities become almost indistinguishable. Thus, 
the perturbative matching between the two smoothers r ~ ric /(3 — 15ci) is confirmed as far as results on the 
topological susceptibility are concerned. 

In Fig. [9] we present the topological susceptibility as a function of the average action density 

defined as the common scale for cooling and the gradient flow. The susceptibility x has been extracted for 
the clover and improved definitions of the topological charge density and computed using the ensembles with 
/? = 1.95 and /3 = 2.10 for our three actions. We observe that for all three actions and for both definitions of 
the topological charge density as well as for ric > 2 we obtain very good agreement. For our most accurate 
calculation using the (3 = 1.95 ensemble, results on y obtained using cooling and gradient flow are in excellent 
agreement, but differ for the clover and improved definitions of x- Complementarily, for our finest lattice 
spacing ensemble with /3 = 2.10, we observe that the topological susceptibilities for the clover and improved 
definitions of the topological charge density become closer for ric — 6, 10, 20 for Wilson, Symanzik tree-level 
improved and Iwasaki smoothing actions, respectively. This is in accordance with the fact that the topological 
susceptibility based on the two different definitions of the topological charge density is expected to become 
the same towards the continuum limit. 

Returning to Fig. [5] one can see that there is a plateau for the topological susceptibility as a function of 
the smoothing scale when the clover/plaquette definitions are used for the topological charge density which 
sets in when the Wilson action is used for ric ^ 40. A plateau is also observed for the improved definition if 
the Symanzik tree-level improved action is used for ric 40. On the contrary, when using the Iwasaki action, 
the susceptibility increases with ric (or equivalently with r). This means that ric is not large enough for the 
Iwasaki action. 
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Figure 9: The susceptibility as a function of the average action density {So) for /3 = 1.95 (top) and /3 = 2.10 
(bottom) ensembles and the three actions. 


5.4 Correlation Coefficient 

In the previous sections we showed that cooling and gradient flow provide results, which are equivalent for the 
average action density and the topological susceptibility under the perturbative rescaling of Eq. (l40l) . In this 
section, we examine the linear correlation coefficient for these two procedures, defined as 


^Ql(nc),Q2(r) 


{{Qi-Qi) {Q2-Q2)) 

y((Qi-Qi)^)((Q2-Q2)^) 


(46) 


where {Qi{nc)} and {Q 2 {nc)} are the two sets of values of the topological charge obtained using cooling at 
ric and gradient flow at r respectively on the same gauge configurations. This implies that is a 

matrix of size Uc x r. The question we would like to answer in this section is the level of correlation between 
sets of topological charges extracted via cooling and gradient flow using the same action. For this discussion 
we employ the topological charge using the improved definition. The results for the other two definitions are 
similar. In Fig. (TU] we represent results for the correlation coefficient using the three actions for our three 
ensembles. 

We show the diagonal elements (for ric = (3 —I5ci)r) of the correlation coefficient matrix when 

ric and r are matched with the perturbative expression Eq. dMl). When the Wilson action is used, we observe 
that for ric > I and as we increase ric the coefficient drops till it reaches a nearly stable value (ric > 10 — 20). 
This value is approximately ^93.5% for /3 = 1.90, ~95% for jS = 1.95 and ^98% for /3 = 2.10. Clearly, as 
a —/ 0 the correlation coefBent approaches unity. This indicates that the correspondence between cooling and 
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Figure 10: The diagonal elements of the correlation coefficient CQ^(^nc),Q 2 {T) defined in Eq. (H51) for the Wilson, 
Symanzik tree-level improved and Iwasaki actions respectively for f3 = 1.90, P = 1.95 and (3 = 2.10. We 
consider the topological charge extracted using the improved definition. 


gradient flow has a well-defined continuum limit. A similar behavior is observed when the Symanzik tree-level 
improved action is used obtaining ~93% for /? = 1.90, ^95% for /3 = 1.95 and ~97.5% for /3 = 2.10. Finally 
and likewise when the Iwasaki smoothing action is used the level of correlation is ~92.5% for /3 = 1.90, ~94% 
for /3 = 1.95 and ~96.5% for (3 = 2.10. 
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Figure 11: The correlation coefficient matrices for the Wilson and Symanzik tree-level improved 

action. We consider topological charge extracted for /3 = 1.95 and the improved definition of the topological 
charge density. 


In Fig. [TT]we provide density plots for the full correlation coefficient matrix CQ^(^nc)Q 2 ('r) foi' the Wilson 
and Symanzik tree-level improved actions obtained when the improved topological charge is employed for 
(3 = 1.95. When excluding the very first cooling steps (e.g. ric < 10), the matrix cq^q^ appears to be nearly 
diagonal with the diagonal line denoting the equation ric = {3 — 15ci)t. This behavior is more pronounced for 
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the case of the Wilson and Symanzik tree-level improved actions. Thus these results corroborate the fact that 
cooling the gauge configurations with ric steps has almost the same effect as evolving these configurations via 
gradient flow for r = nc/(3 — 15ci). We expect that at the continuum limit the corresponding distributions 
become perfectly diagonal with the maximum along the diagonal and corresponding to a correlation coefficient 
of 100%. 

6 Conclusions 

In this article we provide a comparison of the results on observables such as the topological charge and the 
susceptibility obtained using gradient flow or cooling. It extends the analysis of Ref. [7] to include gauge 
actions with rectangular terms. The comparison is realized both analytically in perturbation theory and 
numerically. For our analytic analysis we follow the perturbative treatment of Ref. [7] , which was performed 
for the Wilson action and we show how to generalize it to Symanzik improved actions with rectangular parts. 
More specifically, we derive the corresponding relation between the continuous gradient flow time r and the 
number of the discrete cooling steps ric by expanding the flow steps perturbatively including terms up to 0{a). 
The relation we obtain is r ~ ndii — 15ci) where Ci is the coefficient which multiplies the rectangular term 
in the gauge action. This becomes exact as a —>■ 0 and does not depend on the details of the gauge group; 
although this is derived for SU{2) the generalization to SU{3) is straight forward. For the numerical results 
we use configurations produced with Nf = 2 -|- 1 -|- 1 twisted mass fermions and the Iwasaki gauge action. 
Although strictly speaking the relation we derived is valid only as a —>■ 0, we confirm numerically that the 
action density, used as a common scale, coincides for both procedures. 

By investigating the time histories of the topological charge we observe that these behave in the same 
manner for both smoothing procedures indicating equivalence between them. The histograms of the topo¬ 
logical charge distributions for fixed ric and r ^ rid{3 ~ 15^c) are almost the same for both smoothers and 
approximately Gaussian having the same width. This already suggests an equivalence for the topological sus¬ 
ceptibility, which is confirmed by calculating the topological susceptibility x for all three lattice definitions of 
the topological charge density as a function of the smoothing scale and the average action for both smoothers. 
This enables us to demonstrate that after a very few cooling steps ric ^ 2 the topological susceptibility for gra¬ 
dient flow and cooling become equivalent; this holds for all tested smoothing actions and all lattice definitions 
of the topological charge. 

Finally we look at the correlation coefficient, which can be used to reveal similarities between the different 
definitions of the topological charge. We observe maximum correlation for gauge configurations that have 
been smoothed via gradient flow or cooling according to the relation r ~ rid {3 ~ 15ci). In addition, we 
observe that after a few cooling steps the correlation coefficient becomes stable with increasing value towards 
the unity as we approach the continuum limit (decreasing the lattice spacing). For instance already for our 
finest lattice with [3 = 2.10 the correlation coefficient is ^ 98%, ~ 97.5% and ^ 96.5% when smoothing with 
Wilson, Symanzik tree-level improved and Iwasaki action, respectively. 

The main conclusion of this study is that one can use cooling or the gradient flow in order to extract the 
topological properties of configurations smoothed with gauge actions, which include square and rectangular 
terms. This equivalence is manifested by using the relation r = nd {3 — 15ci) derived in perturbation theory. 
In practice, this means that one may opt to use cooling to extract the topological charge Q. An approximate 
comparison between the gradient flow time r with integration step e = 0.01 and cooling step ric for an action, 
which includes rectangular terms gives cpu_time(r = l)/cpu_time(ric = 1) — 160. Hence, for the Symanzik 
tree-level improved action, gradient flow is slower than cooling by a factor of ~ 38 while for the Iwasaki action 
by ^ 20. These estimates depend on the integrator used for the gradient flow and the integration step e. The 
speed-up cooling gives in comparison to gradient flow is 0(10) and this could decrease the computational cost 
by the same factor in investigations where one is mainly interested in the topological susceptibility and where 
a large number of configurations is required. 
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